We make a self-propelling Purcell Swimmer and dip it in a highly viscous fluid.
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import seaborn as sns; sns.set(); sns.set_style("whitegrid")
plt.rcParams["animation.html"] = "jshtml"
from celluloid import Camera
import stokeslets as slts
import purcell as pur
import matplotlib.patches as ptc
figfile = "figs/"
Purcell imagines a swimmer made of three 'links' or rods connected by two joints. The geometry of the swimmer can be described at any given point in time by the two angles $\theta_1$ and $\theta_2$. It undergoes a cycle as shown below from $S_1$ to $S_5$ moving its arms alternately.

We need to define the position of the rods and the motors. Let's construct Purcell's and place the end points of the rod in the $S_1$ configuration. The rod is parametrised by $\eta$ — the ratio of the body length to the arm length — and $\phi$ = $\max_t(\theta_1(t)) = \max_t(\theta_2(t))$.
et = 2 #ratio of the body to one arm
ph = np.pi/3
size = 1 #total length
s = np.zeros([4,3]);
a = size/(2+et) #one-arm length
b = a*et #body length
s[1,0:2] = np.array([-b/2,0]) #motor 1
s[2,0:2] = np.array([b/2,0]) #motor 2
s[0,0:2] = np.array([-a*np.cos(ph)-b/2,a*np.sin(ph)]) #arm 1
s[3,0:2] = np.array([a*np.cos(ph)+b/2,a*np.sin(ph)]) #arm 2
r = np.array([1,2]) #motor specifications
fig,ax = plt.subplots(figsize=(4,4))
ax.plot(s.T[0],s.T[1],'k-',alpha=0.5)
arc1 = ptc.Arc((s.T[0,2],s.T[1,2]),0.2,0.2,angle=0,
theta1=0,theta2=ph*180/np.pi,color='r',linewidth=1.5)
ax.add_patch(arc1)
plt.annotate(r"$\phi$",(s.T[0,2]+0.15,s.T[1,2]+0.075))
plt.xlim(-1,1); plt.ylim(-1,1);
plt.savefig(figfile+"purcell_swimmer_configuration.png",dpi=300,bbox_inches='tight')
plt.show()
Define the time vector.
dt = 0.01
T = 16
tim = np.arange(0,T,dt)
Let's define $\tau_1 = [0,0,q_1(t)]$ and $\tau_2 = [0,0,q_2(t)]$.
tau = np.array([np.zeros(np.shape(tim)+(3,))]*np.size(r))
In Purcell's paper, the swimmer moves one arm at a given time. This can be achieved by a constant torque difference between the two motors. To describe this motion, we need a parameter $P = |q_2(0)|$ which is the maximum torque exerted by one of the motors. For simplicity, we assume that the other motor exerts no torque when the first is switched on.
P = 1.463 #maximum power
tau = pur.stroke_input(P,tim,r,strokes=2) #the second motor exerts no torque when the first is moving
plt.plot(tim,tau[0,:,-1],'-',alpha=0.9);
plt.plot(tim,tau[1,:,-1],'-',alpha=0.9);
plt.legend(["$q_1 (t)$","$q_2 (t)$"],loc="upper right");
plt.xlabel(r"$t \ (a.u)$"); plt.ylabel(r"$\tau_z \ (a.u)$")
plt.savefig(figfile+"torque_over_time.png",dpi=300,bbox_inches="tight")
R = slts.mesher()
_,pswimmer = slts.evolve(tau,tim,R,r,s,k=100,e=0.3,c=0.6)
fig = plt.figure(figsize = (4.5,4.5))
plt.xlim(-1,1); plt.ylim(-1,1)
camera = Camera(fig)
for i in range(0,np.shape(pswimmer)[0],8):
plt.plot(pswimmer[i].T[0],pswimmer[i].T[1],'k-', linewidth=2, alpha=0.5, label='_nolegend_')
plt.plot(np.mean(pswimmer[i].T[0]),np.mean(pswimmer[i].T[1]),'ro', linewidth=2, alpha=0.5)
camera.snap()
plt.legend(["CoM"])
plt.close()
animation = camera.animate()
animation.save(figfile+"purcell_swimmer.gif",writer = 'imagemagick',fps = 20,dpi=300)
animation
We can try finding the torque difference required for a Purcell stroke by trial and error, but it would be easier to sample different $P$ values and see which corresponds to a Purcell stroke.
Defining $\theta_1$ and $\theta_2$ by Purcell's convention as in the figure above, we can come up with a function that must be equal to zero for a Purcell stroke after half a cycle:
$$\Delta \phi = (\phi-\theta_2)+(\phi+\theta_1)$$print("P0 = ",pur.delphihalf(tau,pswimmer,ph))
P0 = 0.04077787615948769
Now, we can try sampling over the parameter space to see which values of $P$ give $\Delta\phi_{1/2} = 0$ where $\Delta \phi _{1/2}$ is $\Delta \phi$ after half a cycle is completed.
et = 2 #eta value
ph = np.pi/3 #phi value
res = 0.1
P_vec = np.arange(0.6,1.6+res,res) #range of P values
delphiH = np.zeros(np.shape(P_vec))
for i,P in enumerate(P_vec):
s,r,tim = pur.init(et,ph,T=9)
tau = pur.stroke_input(P,tim,r,strokes=9/8)
R = slts.mesher()
_,pswimmer = slts.evolve(tau,tim,R,r,s,k=100,e=0.3,c=0.6)
delphiH[i] = pur.delphihalf(tau,pswimmer,ph)
print("Progress:",np.round((i+1)/np.size(P_vec)*100,2),"%")
Progress: 9.09 % Progress: 18.18 % Progress: 27.27 % Progress: 36.36 % Progress: 45.45 % Progress: 54.55 % Progress: 63.64 % Progress: 72.73 % Progress: 81.82 % Progress: 90.91 % Progress: 100.0 %
plt.figure(figsize=(5,3.5))
plt.plot(P_vec,delphiH,'o',color='seagreen');
plt.ylabel(r"$\Delta\phi_{1/2}$")
plt.xlabel(r"$P$")
plt.savefig(figfile+"delphiH_vs_P_example.png",dpi=300,bbox_inches="tight")
plt.show()
The nice one-to-one mapping simplifies our job by quite a bit. We can use bisection method to find $P_0$.